/* ************************************************************************* **
**    OpenVFIFE - Open System for Vector Form Instrinsic                     **
**                Finite Element Method (VFIFE)                              **
**                GinkGo(Tan Biao)                                           **
**                                                                           **
**                                                                           **
** (C) Copyright 2021, The GinkGo(Tan Biao). All Rights Reserved.            **
**                                                                           **
** Commercial use of this program without express permission of              **
** GinkGo(Tan Biao), is strictly prohibited.  See                            **
** file 'COPYRIGHT'  in main directory for information on usage and          **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.                  **
**                                                                           **
** Developed by:                                                             **
**      Tan Biao (ginkgoltd@outlook.com)                                     **
**                                                                           **
** ************************************************************************* */

// $Date: 2021-06-07$
// Written: Tan Biao
// Revised:
//
// Purpose: This file contains the class definition for elements, including
// base class BaseElement (abstract base class), StructElement, PlaneElement,
// SolidElement, and several elements derived from StructElement - Link2D,
// Link2DLD, Link3D, Link3DLD, Cable3D, Beam2D and Beam3D.

// The interface:
//

#ifndef ELEMENT_H_
#define ELEMENT_H_
#include <vector>
#include <map>
// #include <vector>
#include <fstream>
#include "Eigen/Dense"
#include "particle.h"
#include "material.h"
#include "section.h"
#include "gridsection.h"
#include "integrator.h"

typedef std::array<double, 3> StdArray3d;
struct Strain
{
    /* organize strain at one point, for outputing */
    double sx, sy, sz;
    double sxy, syz, sxz;
};

struct Stress
{
    /* organize stress at one point, for outputing */
    double sx, sy, sz;
    double sxy, syz, sxz;
};

struct ElementForce
{
    /* organize force of rod element, for outputing */
    double fx;
    double sfy, sfz;
    double mx, my, mz;
};

/* Sorry, I have no time to write a detailed comments for every element class.
* In fact, I am not very satisfied with the current framework. I will redesign
* openVFIFE after graduating. When the time comes, a detailed documentation
* will be released with new version of openVFIFE. Sorry again, there will be
* no comments of element classes.
*
*/
class BaseElement
{
    protected:
        int id_;
        std::string type_;
        int num_of_particles_;
        std::map<int, Particle*> particles_;
        const BaseMaterial* material_;
        const BaseSection* section_;
        double mass_;
        Strain strain_;
        Stress stress_;
        Eigen::Vector3d ex_, ey_, ez_;
        bool is_failed_;

    protected:
        inline bool checkParticles();
        virtual void setParticleDof()=0;
        virtual void setParticleMass()=0;
        virtual void calcOrientVector()=0;
        virtual void calcElementForce()=0;

    public:
        BaseElement(int id);
        virtual ~BaseElement();
        void addParticle(Particle* p);
        void setMaterial(BaseMaterial* mat);
        void setSection(BaseSection* sect);

        virtual void info() const;
        int id() const;
        std::string type() const;
        int num_of_particles() const;
        const std::map<int, Particle*>* particles() const;
        const BaseMaterial* material() const;
        const BaseSection* section() const;
        double mass() const;
        const Eigen::Vector3d* ex() const;
        const Eigen::Vector3d* ey() const;
        const Eigen::Vector3d* ez() const;
        const Strain* strain() const;
        const Stress* stress() const;
        bool is_failed() const;
        virtual double calcCriticalTimeStep()=0;
        virtual void setParticleForce()=0;
        virtual void outputForce(std::fstream &fout)=0;
        virtual void outputStress(std::fstream &fout)=0;
        virtual void outputStrain(std::fstream &fout)=0;
};


class StructElement: public BaseElement
{
    private:
        inline void calcLength();
        inline void calcMass();

    protected:
        Particle* pi_;
        Particle* pj_;
        StdArray6d fi_;
        StdArray6d fj_;
        double length_;
        double lt_, la_;
        ElementForce element_force_;

        virtual void setParticleDof()=0;
        virtual void setParticleMass()=0;
        virtual void calcOrientVector()=0;
        virtual void calcElementForce()=0;

    public:
        StructElement(int id, Particle* pi, Particle* pj,
            const BaseMaterial* mat, const BaseSection* sect);
        virtual ~StructElement();
        double length() const;
        virtual double calcCriticalTimeStep();
        virtual void setParticleForce();
        // virtual void kill()=0;
        virtual void outputForce(std::fstream &fout);
        virtual void outputStress(std::fstream &fout);
        virtual void outputStrain(std::fstream &fout);
};


class PlaneElement: public BaseElement
{
    // TODO:
    private:

    protected:
        virtual void setParticleDof()=0;
        virtual void setParticleMass()=0;
        virtual void calcOrientVector()=0;
        virtual void calcElementForce()=0;

    public:
        PlaneElement(int id, std::vector<Particle*> ps, BaseMaterial* mat,
            BaseSection* sect);
        virtual ~PlaneElement();
        virtual double calcCriticalTimeStep()=0;
        virtual void setParticleForce()=0;
        virtual void outputForce(std::fstream &fout)=0;
        virtual void outputStress(std::fstream &fout)=0;
        virtual void outputStrain(std::fstream &fout)=0;
};


class SolidElement: public BaseElement
{
    // TODO:
    private:

    protected:
        virtual void setParticleDof()=0;
        virtual void setParticleMass()=0;
        virtual void calcOrientVector()=0;
        virtual void calcElementForce()=0;

    public:
        SolidElement(int id, std::vector<Particle*> ps, const BaseMaterial* mat,
            const BaseSection* sect);
        virtual ~SolidElement();
        virtual double calcCriticalTimeStep()=0;
        virtual void setParticleForce()=0;
        virtual void outputForce(std::fstream &fout)=0;
        virtual void outputStress(std::fstream &fout)=0;
        virtual void outputStrain(std::fstream &fout)=0;
};


class Link2D:public StructElement
{
    protected:
        const StdArray6d* posi_;
        const StdArray6d* posj_;

        virtual void setParticleDof();
        virtual void setParticleMass();
        virtual void calcOrientVector();
        virtual void calcElementForce();

    public:
        Link2D(int id, Particle* pi, Particle* pj, const BaseMaterial* mat,
               const BaseSection* sect);
        virtual ~Link2D();
        virtual void outputForce(std::fstream &fout);
};


class Link2DLD:public Link2D
{
    private:
        /* parameters for plasitc material
        * ds - delta strain;
        * stress - current stress,
        * strain - current strain,
        * ps - current plastic strain
        * alpha - back stress,
        * q - hardening funciton,
        */
        double mstress_, mstrain_, malpha_, mps_, mq_;

    protected:
        virtual void calcOrientVector();
        virtual void calcElementForce();

    public:
        Link2DLD(int id, Particle* pi, Particle* pj, const BaseMaterial* mat,
                const BaseSection* sect);
        virtual ~Link2DLD();
};


class Link3D: public StructElement
{
    protected:
        const StdArray6d* posi_;
        const StdArray6d* posj_;

        virtual void setParticleDof();
        virtual void setParticleMass();
        virtual void calcOrientVector();
        virtual void calcElementForce();

    public:
        Link3D(int id, Particle *pi, Particle *pj, const BaseMaterial *mat,
            const BaseSection *sect);
        virtual ~Link3D();
        virtual void outputForce(std::fstream &fout);
};


class Link3DLD:public Link3D
{
    private:
        /* parameters for plasitc material
        * ds - delta strain;
        * stress - current stress,
        * strain - current strain,
        * ps - current plastic strain
        * alpha - back stress,
        * q - hardening funciton,
        */
        double mstress_, mstrain_, malpha_, mps_, mq_;

    protected:
        virtual void calcOrientVector();
        virtual void calcElementForce();

    public:
        Link3DLD(int id, Particle* pi, Particle* pj, const BaseMaterial* mat,
               const BaseSection* sect);
        virtual ~Link3DLD();
};


class Cable3D:public Link3DLD
{
    /* Space Cable element with prestress, only support Elastic material */
    protected:
        double prestrain_;
        virtual void calcElementForce();

    public:
        Cable3D(int id, Particle* pi, Particle* pj, const BaseMaterial* mat,
                const BaseSection* sect, double prestrain=0.0);
        virtual ~Cable3D();
};

class Beam2D:public StructElement
{
    /* Elastic 2D Beam element, which can count on geometric nonlinearity
    * Only Elastic materials are supported, if an plasticity material is given
    * only the elastic state will be used. Obviously this element do not
    * support the failure of element.
    */
    protected:
        const StdArray6d* curr_posi_;
        const StdArray6d* curr_posj_;
        const StdArray6d* prev_posi_;
        const StdArray6d* prev_posj_;

        Eigen::Matrix3d stiffness_;
        Eigen::Vector3d force_;

        ElementForce emfi_, emfj_;
        Eigen::Vector3d exa_, eya_, eza_;
        double elem_rotation_angle_;

        inline void initiateStiffnessMatrix();
        virtual void setParticleDof();
        virtual void setParticleMass();
        virtual void calcOrientVector();
        virtual void calcElementForce();

    public:
        Beam2D(int id, Particle* pi, Particle* pj, const BaseMaterial* mat,
               const BaseSection* sect);
        virtual ~Beam2D();
        virtual double calcCriticalTimeStep();
        virtual void outputForce(std::fstream &fout);
        virtual void outputStress(std::fstream &fout);
        virtual void outputStrain(std::fstream &fout);
};

class Beam3D:public StructElement
{
    /* Elastic 3D Beam element, which can count on geometric nonlinearity
    * Only Elastic materials are supported, if an plasticity material is given
    * only the elastic state will be used. Obviously this element do not
    * support the failure of element.
    */
    protected:
        StdArray6d reference_node_;
        const StdArray6d* curr_posi_;
        const StdArray6d* curr_posj_;
        const StdArray6d* prev_posi_;
        const StdArray6d* prev_posj_;

        Eigen::Vector3d exa_, eya_, eza_;
        Eigen::Vector3d gamma_;
        Eigen::Matrix3d r_gamma_;
        Eigen::Matrix3d omega_;
        Eigen::Vector3d em1h_, em2h_, ef1h_, ef2h_;
        ElementForce emfi_, emfj_;

        inline void initiateBeam3D();
        void initiatePrincipleCoord();
        void rotatePrincipleAxial();

        virtual void setParticleDof();
        virtual void setParticleMass();
        virtual void calcOrientVector();
        virtual void calcElementForce();

    public:
        // create spatial Beam with default reference node (0, 0, 1)
        Beam3D(int id, Particle* pi, Particle* pj, const BaseMaterial* mat,
               const BaseSection* sect);
        // create spatial Beam with giving reference node
        Beam3D(int id, Particle* pi, Particle* pj, const BaseMaterial* mat,
               const BaseSection* sect, const StdArray3d ref_node);
        virtual ~Beam3D();
        virtual double calcCriticalTimeStep();
        virtual void outputForce(std::fstream &fout);
        virtual void outputStress(std::fstream &fout);
        virtual void outputStrain(std::fstream &fout);
};


class PlasticBeam2D:public Beam2D
{
    private:
        int integrals_points_;
        GaussLobatto lobatto;
        std::vector<SectionResponse*> response_;
        const SectionGrid* grid_;

    protected:
        inline void isFailed();
        virtual void calcElementForce();

    public:
        PlasticBeam2D(int id, Particle* pi, Particle* pj,
                     const BaseMaterial* mat, const BaseSection* sect,
                     const SectionGridParameter &para);
        virtual ~PlasticBeam2D();
};


class PlasticBeam3D: public Beam3D
{
    private:
        int integrals_points_;
        GaussLobatto lobatto;
        std::vector<SectionResponse*> response_;
        const SectionGrid* grid_;

        inline void initiatePlasticBeam3D(const BaseMaterial* mat,
            const BaseSection* sect, const SectionGridParameter &p);

    protected:
        inline void isFailed();
        virtual void calcElementForce();

    public:
        // create spatial Beam with default reference node (0, 0, 1)
        PlasticBeam3D(int id, Particle* pi, Particle* pj,
            const BaseMaterial* mat, const BaseSection* sect,
            const SectionGridParameter &para);
        // create spatial Beam with giving reference node
        PlasticBeam3D(int id, Particle* pi, Particle* pj,
            const BaseMaterial* mat, const BaseSection* sect,
            const StdArray3d ref_node, const SectionGridParameter &para);
        virtual ~PlasticBeam3D();
};


#endif